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Abstract 

The theoretical dispersion curves n{uj) (refractive index versus fre¬ 
quency) of ionic crystals in the infrared domain are expressed, within 
the Green-Kubo theory, in terms of a time correlation function in¬ 
volving the motion of the ions only. The aim of this paper is to inves¬ 
tigate how well the experimental data are reproduced by a classical 
approximation of the theory, in which the time correlation functions 
are expressed in terms of the ions orbits. We report the results of 
molecular dynamics (MD) simulations for the ions motions of a LiF 
lattice of 4096 ions at room temperature. The theoretical curves thus 
obtained are in surprisingly good agreement with the experimental 
data, essentially over the whole infrared domain. This shows that 
at room temperature the motion of the ions develops essentially in a 
classical regime. 
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Experimental data on the dispersion curves n{u) - refractive index ver¬ 
sus angular frequency - of ionic crystals in the infrared domain have been 
available for a long time [1]. However, a microscopic explanation covering 
the whole infrared domain is apparently lacking, notwithstanding the fact 
that the general theoretical frame is well dehned. On the other hand, a re¬ 
cent work [2] (see also [3]) has shown that the dispersion relation a;(k) in 
ferroelectrics can be calculated through MD simulations. In this paper we 
show that a calculation of the ions motions through MD simulations leads to 
theoretical dispersion curves n{uj) for ionic crystals which agree in a surpris¬ 
ingly good way with the experimental data essentially over the whole infrared 
domain. In fact, the problem can be dealt with through MD simulations be¬ 
cause by Green-Kubo theory the refractive index is expressed in terms of a 
time correlation function involving the microscopic dipole moment P due to 
the ions motions. 

Let us recall the theoretical frame. The refractive index is just the square 
root of the electric permittivity tensor e = eij{u), which in turn is related to 
the dielectric susceptibility tensor x{^) through 


eij{uj) = 6ij + Anxijiuj) ■ 


( 1 ) 


Now, the susceptibility tensor x{^) is the response function of the considered 
medium to an external electric held, so that the Green-Kubo theory provides 
for it a formula involving a microscopic quantity, the polarization P(t). For 
a ionic crystal the contribution to polarization in the infrared range comes 
mainly from the ions. So the microscopic polarization can be written as 



( 2 ) 


where qs,h is the displacement of the ion of the s species (having charge e^) in 
the h-th cell from its equilibrium position, and the summation is understood 
over all ions belonging to the (small) volume V. 

The quantum theoretical formula for the ionic contribution to suscepti¬ 
bility is given for instance in |1] (see formula (2.16a)) and in [5]. For a review 
see [6]. The corresponding classical limit is then (see formula (2.16b) in [1], 
or see [7| for a purely classical deduction) 



( 3 ) 
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Here, as usual, f3 = l/ksT is inverse temperature with ks the Boltzmann 
constant, while (...) denotes ensemble average. Instead, the contribution of 
the electrons to susceptibility in the infrared is well known from experiments 
to just reduce to a constant (isotropic) tensor which can be taken 

from pp. 

Thus the computation of the refractive index in the infrared range is 
reduced to a computation of the ionic susceptibility, which in turn is reduced, 
through ^ and ([^, to a determination of the classical motions of the ions. 

So, we model the crystal as a face centered cubic lattice of 2N particles 
(the ions, with their atomic masses) in a cubic box of side L = /o(iV/4)3, 
where /q is the lattice constant^ The ions interact through mutual electric 
forces, and through a phenomenological short-range repulsive pair potential 
Vrep- Periodic boundary conditions are imposed, and the equilibrium posi¬ 
tions of the atomic species s = 1, 2 in the cell h G of the lattice, as 
well as the value of the lattice constant Iq, are taken from the literature [8]. 
In the numerical simulations, we took 2N = 4096. 

For the aim of estimating the electric forces between ions, it is known 
(see [9]) that the ions can be dealt with as point particles only if a suitable 
“effective charge” Ce// is substituted for the ions charge. In this paper we 
concentrate on lithium fluoride, which is the halide presenting the simplest 
infrared spectrum, using for it the value given in the literature (see m), 
namely e^ff = 0.81e, where e is the electron charge. 

Concerning the phenomenological repulsive short-range pair potential 
Vrep, in the literature one often takes it to be spherically symmetric, indepen¬ 
dent of the ions’ species, and depending on the distance r as an exponential, 
aexp(—r/r*). Here we chose the analytic form previously considered by Born 
in his MIT lectures m and in the paper [T2|, namely, Vrep{r) = ar p, with 
a cut-off at 5 A. Thus, having hxed the effective charge, the model still con¬ 
tains two free parameters, a and p. We are conhdent that the choice of the 
repulsive potential, exponential or inverse power, doesn’t play an essential 
role. 

The two parameters a and p can in principle be determined through the 
thermodynamic available data, for example from the compressibility and the 
heat of formation. Here we determined them from the experimental disper¬ 
sion relation a;(k) obtained from neutron scattering. To this end we made 

^Recall that the primitive cell of a fee crystal in not cubic at all but tetrahedral, and 
that the unit cubic cell is arranged using four primitive cells. 
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Figure 1: Real part of the refractive index versus frequency for a LiF crystal 
at 300 K. Circles and diamonds are experimental data taken from [T] and [15] 
respectively. Solid line is obtained, through Green-Kubo theory, from MD 
simulations for the motions of a system of 4096 ions. The inset is a zoom of 
the leftmost part of the hgure. 
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reference to the work [12], which was devoted to an analytical estimate of 
the dispersion relation in the linear approximation with the aim of a micro¬ 
scopic explanation of the existence of polaritons. In that work the elastic 
constants of the crystal had to be determined by £t with the experimental 
curves. As the elastic constants are related to the first two derivatives of the 
phenomenological potential V^ep (see [I3|), this information, together with 
the choice that p be an integer, allowed to determine a = 1.25 and p = 6. 

Concerning the electric forces, due to their long-range nature a prob¬ 
lem arises in writing down explicitly the equations of motion when periodic 
boundary conditions are imposed, because one then has to determine the 
electric held produced by an inhnite lattice of charges. As is well known 
(see lEl). the problem is solved in a computationally efficient way using the 
Ewald resummation formula. This amounts to introducing the potential 
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where erfc(a;) is the complementary error function, 5 is a parameter to be 
chosen to insure fast convergences of the series, = (xs^h + qs,h) — (x^/^h' + 
qs',h') and the primed sum over n means that the term for n = 0 is omitted 
if s = s' and h = h'. We took 6 = 0.2 which allows to truncate the 
series at \n^\, \ny\, \n^\ < N^ax = 1 and \k^\, \ky\, < K^ax = 3. 

The numerical integrations of the equations of motion were performed 
with a standard symplectic Verlet method. The integration step was taken 
equal to 2 fs, the duration of each simulation was of 100 ps, and the time 
averages (which will be mentioned in a moment) were taken over 50 ps. 
The initial data for each orbit were assigned by setting the particles in their 
equilibrium positions, while the velocities were extracted according to the 
Maxwell-Boltzmann distribution, with the constraint that the center of mass 
of the system has vanishing velocity. Then, the temperature of the system 
(300 K) was determined through the mean kinetic energy, once thermaliza- 
tion had been attained. 

So we computed the orbits, i.e. the displacement qs,h(^) of the ions from 
their equilibrium positions, for a LiF lattice of 4096 ions at 300 K, and thus 
the microscopic polarization P(t) too, given by (|^. In formula ([^ for the 
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Figure 2: Imaginary part of the refractive index versus frequency for a 
LiF crystal at 300 K. As in figure [TJ circles are experimental data taken 
from [1], while solid line is obtained, through Green-Kubo theory, from MD 
simulations for the motions of a system of 4096 ions. 
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Figure 3: Imaginary part of the electric permittivity £ (in logarithmic scale) 
versus rescaled frequency X = oj/col, where col = 665 cm“^ is the frequency 
of the zero longitudinal optical mode. Circles are experimental data from [T], 
while solid line is obtained, through Green-Kubo theory, from MD simula¬ 
tions for the motions of a system of 4096 ions. Compare with the theoretical 
curve reported in Fig. 12 of ref. lai. obtained with quantum mechanical 
methods. Note that the experimental data reported there are the older ones 
[22] and [23] • The agreement with the new data [T] reported here is appar¬ 
ently better. 
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susceptibility the ensemble average was estimated as the mean value, over 10 
different orbits, of the corresponding time average. 

As a check for our computations we controlled that the susceptibility ten¬ 
sor actually is, as expected for LiF, an isotropic one. Namely, the off-diagonal 
matrix elements are negligible and the diagonal ones can be considered as 
equal. One may presume that these properties should not occur in a less 
symmetric type of crystal, so that phenomena of optical activity would show 
up in such cases. In our case, however, the susceptibility tensor can be dealt 
with as a scalar, which was estimated as the mean of the three diagonal 
matrix elements. 

The results of the computations are now illustrated. In figure we show 
the main result of this paper, namely, the real part of the refractive index n = 
^/e, with e given by Q and x computed through (|^. The theoretical curve 
was plotted together with the experimental data taken from [1] and from 
|15j . The error bars of the experimental data are not reported because the 
relative errors are negligible (less than 0.1 %). Notice that, in the theoretical 
curve, only the parameter was directly fitted to the set of experimental 
data, taking it from the data reported in [1]. 

The overall agreement appears to be surprisingly good, essentially over 
the whole infrared domain. Actually one may notice that, as exhibited by 
the inset, apparently the agreement is not that good in the far infrared range, 
particularly below 100 cm“^. In this connection we point out first of all that 
in such a region the experimental data appear not to be so sure because in 
the literature different series of measures are found, which are not mutually 
compatible in view of the declared errors. Our theoretical curve are however 
in fairly good agreement with the experimental data of reference [I^ , which 
are, at any rate, the more recent ones. 

Should the measurements of ref. na be the good ones, then a very in¬ 
teresting perspective would be opened. Indeed it is known that for very low 
frequencies (below 1 MHz) the refractive index is slightly above 3 (between 
3.05 and 3.018 according to [IS]), while in the far infrared range we found, in 
agreement with ref. [IS], a value near 2.9, i.e. below the static limit a; —)■ 0. 
So, the refractive index should increase as u tends to zero, and this indeed 
happens in the microwave range (GHz), according to [H], where a value 
slightly below 3 is found. 

To settle this question, at least at the level of MD simulations, integrations 
of the equations of motion over much longer time scales would be necessary. 
By the way, the problem of the dependence of the results on the Integra- 


tion time is a very delicate one, and was discussed, for example, in [TS] in 
connection with a MD computation of the specihc heat of an Argon crystal, 
where a comparison between MD simulations and experimental data was per¬ 
formed somehow in the same spirit of the present paper. Another possibility 
to tackle the problem is to enlarge the model, in order to take some further 
physical property into account. We are thinking of the retarded character of 
the electromagnetic forces. In fact, in the work 113 retardation turned out to 
be the essential qualitative ingredient in proving the existence of polaritons 
(i.e., the presence of two new branches in the dispersion relation; see [19] and 
pup in a microscopic model. We are analogously suggesting that retardation 
might lead to the arising of new absorption bands in the far infrared range. 
We hope to come back to these problems in the future. 

In any case, leaving aside the problems related to the far infrared, the 
present results show that the classical motions of the ions at room tempera¬ 
ture, estimated by MD simulations on a sample of 4096 ions at 300 K, lead to 
theoretical dispersion curves which are, overall, in very good agreement with 
the experimental data. We recall that theoretical estimates of the dispersion 
curves of ionic crystals were given also in the frame of quantum mechanics 
for the same model. To our knowledge such theoretical estimates, which 
use many-phonon perturbation methods, presently do not cover the whole 
infrared domain. In the domain in which they apply, the results are appar¬ 
ently not better than the present ones. This is illustrated through figure 
which should be compared with the quantum theoretical predictions for LiF 
at room temperaturereported in Fig. 12 of 1211 . 
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